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ABSTRACT 


This  paper  describes  a procedure  for  computing  the  maximum  likelihood 
estimates  of  the  parameters  of  the  distribution  of  the  sum  of  three  inde- 
pendent exponential  random  variables.  By  fitting  sample  interevent  time 
da*a  from  a real  system  to  this  distribution,  one  can  create  a simulation 
of  the  system  that  exploits  the  regenerative  representation  of  queueing 
systems  [3]  to  analyze  the  simulation's  output  by  relatively  elementary 
statistical  methods.  The  paper  also  describes  computation  of  the  sample 
asymptotic  covariance  matrix  and  an  implementation  of  the  likelihood  ratio 
for  testing  six  hypotheses  that  are  special  cases  of  interest.  A set  of 
FORTRAN  subroutines  for  executing  these  procedures  appears  in  the  Appendix. 
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1.  Introduction 

The  purpose  of  this  paper  is  to  describe  a procedure  for  computing 
the  maximum  likelihood  estimates  (MLE)  of  the  parameters  of  the  distri- 
bution of  the  sum  of  three  independent  exponentially  distributed  random 
variables.  Although  the  case  of  equal  parameters  yields  an  Erlang  distri- 
bution, for  which  the  MLE  are  known,  the  more  general  case  has  received 
little  attention  in  the  statistical  literature.  Two  possible  reasons  for 
this  omission  occur  to  the  writer.  Firstly,  since  the  corresponding  MLE 
equations  are  not  amenable  to  analytical  solution,  one  needs  to  employ 
numerical  analytic  techniques  to  solve  them.  Conceptually,  the  presence 
of  multiple  maxima  makes  this  an  onerous  approach.  Secondly,  since  the 
distribution  has  three  parameters,  the  principle  of  parsimony  encourages 
one  to  use  alternative  two  parameter  distributions  whenever  a fit  of  equal 
or  almost  equal  quality  can  be  obtained.  These  distributions  include  the 
garana,  lognormal  and  Weibull.  Choi  and  Wette  [1]  describe  a procedure  for 
computing  the  gamma  MLE.  Thoman,  Bain  and  Antle  [13]  describe  a procedure 
for  computing  the  Weibull  MLE.  Although  both  procedures  rely  on  the 
Newton-Raphson  Iterative  method  no  unusual  problems  arise.  For  the  log- 
normal distribution  the  MLE  relate  directly  to  the  MLE  for  the  corresponding 
normal  distribution.  Johnson  and  Kotz  [9]  discuss  issues  related  to  the 
MLE  for  these  distributions,  including  bias  removal. 

Given  the  attractions  of  alternative  distributions,  a relatively  strong 
justification  for  pursuing  the  research  presented  here  seems  in  order. 

Recent  developments  in  the  field  of  discrete  event  simulation  provide  this 
justification.  In  [5,6]  Fishman  points  out  that  In  the  simulation  of 


queueing  systems  one  could  use  the  exit  of  the  system  from  the  empty  and 
idle  state  to  demarcate  the  sample  path  of  a stochastic  process  of  interest 
into  independent  segments  each  of  which  obeys  the  same  probability  law. 

This  demarcation  enables  one  to  use  relatively  elementary  statistical 
methods  to  compute  point  and  interval  estimates  for  population  parameters 
of  interest  [5,6].  The  most  appealing  theoretical  feature  of  this  obser- 
vation is  that  the  i.i.d.  property  holds  regardless  of  the  distributions 
of  interarrival  and  service  times.  The  most  unappealing  feature  arises 
when  either  the  activity  level  increases  or  the  number  of  servers  increases 
for  a given  activity  level.  In  particular,  the  frequency  with  which  the 
system  exits  the  empty  and  idle  state  declines  dramatically.  In  turn, 
this  can  result  in  excessively  long  simulation  runs  if  one  is  determined 
to  collect  a prespecified  number  of  i.i.d.  segments. 

In  [2]  and  [3]  Crane  and  Iglehart  introduce  the  more  general  notion 
of  a regenerative  process  into  the  analysis  of  simulation  output.  In 
particular,  any  state  can  serve  as  a demarcating  state,  provided  that 
statistical  behavior  after  entry  into  that  state  is  independent  of  behavior 
prior  to  entry  and  that  the  state  occurs  infinitely  often.  States  with 
these  properties  are  called  regenerative.  If  one  can  identify  all  such 
states  then  one  can  use  the  most  frequently  occurring  one  to  demarcate  the 
specified  number  of  i.i.d.  segments.  If  the  Interevent  distributions  are 
exponential  then  all  states  can  serve  this  demarcating  purpose.  Since 
exponential ity  is  too  restrictive  an  assumption  in  general.  Crane  and 
Iglehart  [4]  attempt  to  Identify  approximate  regenerative  states.  Their 
procedure  calls  for  a careful  scrutiny  of  the  particular  system  being 
simulated. 
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* An  alternative  approach  to  realizing  the  regenerative  property  arises 
when  interevent  times  have  continuous  uni modal  distributions.  Then  a 
theoretical  basis  exists  for  approximating  each  of  these  distributions 

by  the  distribution  of  the  sum  of  an  arbitrary  nimber  of  independent 
exponential  random  variables.  In  particular,  one  way  to  look  at  this  is 
to  consider  the  polynomial  approximation  to  the  corresponding  characteristic 
function  where  the  reciprocals  of  the  roots  of  the  polynomial,  which  are  real 
for  unimodality,  are  the  means  of  the  exponential  random  variables. + If  one 
adopts  this  characterization  then  interevent  times  in  the  simulation  become  sums 
of  independent  exponential  random  variables.  Suppose,  interarrival  times  are 

I 

representable  as  the  sum  of  two  independent  exponential  random  variables 

* and  service  times  are  exponential.  Then  by  adding  a new  entry  to  the 
state  vector  that  characterizes  which  of  the  two  stages  the  next  arrival 
occupies,  one  provides  the  mechanism  for  realizing  regenerative  states. 

If  service  times  are  representable  as  the  sum  of  three  independent  ex- 
ponential random  variables  then  three  additional  entries  in  the  state  vector 
to  keep  track  of  the  number  of  jobs  In  each  stage  enable  one  to  exploit  the 
regenerative  property  again.  The  price  paid  for  this  ability  is  the 
Increased  bookkeeping  for  the  state  vector,  an  effic  ?tit  approach  to  which 
is  described  in  [7]. 

Although  the  foregoing  discussion  motivates  the  use  of  distributions 
of  sums  of  independent  exponentials,  a procedure  for 
implementing  the  approach  is  practice  remains  to  be  developed.  Ideally, 
one  would  like  to  fit  such  a distribution  by  the  distribution  of  the  sum 

* of  a large  number  of  exponential  variates  and*'  through  a formal  hypothesis 
testing  procedure,  reduce  that  sum  to  the  minimal  number  necessary  to 

fThis  assumes  a polynomial  in  iu  where  i * /T  . 
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account  for  variation  in  the  data.  The  present  paper  describes  a first 
step  in  this  direction  in  Section  2 by  fitting  the  sum  of  three  i ndependent  ex- 
ponential random  variables  and  then  testing  six  hypotheses  designed  to 
reduce  the  length  of  the  state  vector.  In  particular.  Section  2 describes  a 
procedure  for  finding  the  MLE, their  sample  asymptotic  covariance  matrix 
and  for  using  the  likelihood  ratio  to  test  hypotheses.  The  steps  outlined 
in  Section  2 are  implemented  in  a set  of  FORTRAN  subroutines  in  the 
Appendix. 


-5- 


2 . The  Procedure 


Let  Yp  Yg  and  Y^  be  independent  random  variables  from  E(a), 
E{b)  and  E(c),  respectively,  where  E(e)  denotes  the  exponential 
distribution 


f e-*/8/e 


0)  f(*)  - 


Osxs®  0 < 0 
elsewhere. 


Then  X = Yj  + y2  + Y^  has  the  probability  density  function  (p.d.f.) 
(2)  f(x,a,b,c)  = g(x,a,b,c)  + g(x.b.a.c)  + g(x,c,a,b) 


where 

i'i.  g(x,o,*tp)  « oe^/Co^Ke-p). 

Given  a sample  X^,...,Xn  from  (2),  we  wish  to  compute  §,  b»  c,  the 
MLE  of  a,  b and  c,  respectively.  These  follow  from  maximization  of 
the  likelihood  function 
n 

(4)  L ° n f(X,,a,b,c)  . 

i3l  1 

Here  a,  b,  c asymptotically  have  the  trivariate  normal  distribution 
with  means  a,  b,  and  c,  respectively.and  the  minimum  variance  covariance 
uiatrix  l , where  [10] 
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To  obtain  the  MtE  one  usually  solves 


(6) 


9 In  L _ 
90 


9f(Xi ,a,b,c) 


1=1  f(X^,a,b,c) 


° 0 o a a,b,c 


90 


simultaneously  for  a,  b and  c.  In  the  present  case  (6)  does  not  admit 
an  analytical  solution.  Moreover,  the  only  sufficient  statistics  appear 
to  be  X^,...,Xn  which  do  little  to  ease  the  computational  burden  of  a 
numerical  solution. 


feasible  Region 

Although  the  possibility  of  multiple  maxima  makes  maximization  of  L 
difficult  in  general,  we  can  reduce  some  of  this  difficulty  by  noting 
that 

(7)  f(x,a,b,c)  « f(x,a,c,b)  = f(x,b»a,c) 

= f(x,b,c,a)  = f(x,c.a,b) 
f (x,c,b,a)  • 
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This  Implies  that  L has  at  least  6 maxima  of  equal  magnitude.  Intro- 
ducing the  constraints 

(8)  a s b 5 c 

removes  this  ambiguity.  One  can  also  show  that 

(9)  a2  USLL  + b2  iinj,+  c2  ±t*L  m Q 

aa  ab  sc 

leads  to 

(10)  a + b + c » I 

y a ( Vn)  £ X,  . 

iDl  1 

Now  the  constraints  (8)  and  00}  imply 

(11)  Q s 2a  s 5f  - c 
5f  - a s 2c  s I 

which  define  the  feasible  region  in  the  a -c  space  of  Figure  1 where 
maximization  of  t is  to  occur. 

The  arcs  and  nodes  of  the  feasible  region  in  Figure  1 play  a 
special  role  here.  In  particular,  arcs  AB,  BC  and  AC  correspond  to 
hypotheses  1,  2 and  3 in  Table  1 and  nodes  B,  A and  C correspond  to  the 
Erlang  nypotheses  4,  5 and  6.  In  addition  to  examining  these  special  cases 
in  the  process  of  maximization  of  l,  one  can  use  the  likelihood  ratio  test 
to  evaluate  the  effect  of  assuming  that  one  of  these  special  cases  repre- 
sents the  underlying  structure  of  f in  (1).  This  issue  is  discussed  shortly. 


col  xl 


angle  In  Figure  1.  The  set  of  FORTRAN  subprogram  listed  In  the  Appendix 

* 

effects  a grid  search  on  a in  user  specified  increments  of  6 over 

* a 

tO*  3/3]  and  for  each  a performs  a binary  search  for  c in 


Table  1 
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Distributions  and  Derivatives  Under  Alternative  Hypotheses 


i 

Hi 

f(x,a,b,c) 

af(x,a,b,c) 

a f(x,a0btc) 

a f(x,a,b,c) 

da 

a b 

3 C 

0 

asbsc 

g(x,a,b,c)+g(x,b,a,c)+g(x,c,a,b) 

h-|  (x  ,a  ,b»c) 

■ 

h-|(x,D,a,c) 

; 

h^x.c.a.b) 

1 

a=0,b<;c 

g(x,b,c,0)+g(x,c,b,0) 

- 

(x,b,c,0 ) 

(x,c,b,0) 

B 

a=bsc 

p 

g(x,a,c,c)[x(a-c)-ac]/a  +g(x„c,a,a) 

h3(x,c,a) 

- 

h2(x ,a »c) 

E 

asb=c 

2 

g(x,c,a»a)[x(c-a)-ac]/c  +g(x,a,c,c) 

h2(x,c,a) 

- 

h^fx jS ,c) 

E 

a=ba0<c 

g(x,c,0,0) 

- 

- 

g(x,c,0,0)(x/c-l )/c 

« 

5 

as0,bac 

xg(x,c,0,0)/c 

- 

xg(x,c,0,0)(x/o2)/c2 

4 

6 

x2g(x,c,0,0)/2c2 

- 

1 

x2g(x,e,0,0)(x/c-3)/2c 

hjUtO.G.p^gtx.o.^pKl/o+x/o^l/Ce-sH/io-pJJ+gU^.o.pJ/U-oJ+gfx.p.o.sJ/U-o) 

h2(x,o,p)!3g(x,p,ota)[l/p+x/p^-2/(p-o}]+g(x,o,p,p)[{p-0)x+e(o+p)]/(p-d}8^ 


h,(x,0,o)og(x,p,8,0)(t(o-e)x-8p3[x/p2“l/p-2/(p-o)3+x-0)/pZ+  2g(x,o,p,p)/(o-p) 
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[(X-a)/2,X-2a]  to  within  the  tolerance  6 . The  search  for  c solves 

A 

8£n  L/ac  = 0.  Expression  (10)  yields  b and  the  search  for  the  maximum 
is  effected  by  computation  and  comparison  of  £n  L for  each  computed  set  of 

A A A 

a,  b and  c.  Since  substantial  experience  with  the  UPDATE  subroutine  using 
a complete  grid  search  in  ESTIMA  failed  to  reveal  more  than  one  maximum, 
ESTIMA  was  modified  to  terminate  the  search  once  a maximum  has  been  found. 

The  ARC  and  NODE  subroutines  enable  one  to  check  the  arcs  AB,  BC  and 
AC  and  the  nodes  A,  B and  C for  solutions  that  might  give  improvement. 

Also  HYP123  and  ERLANG  use  the  results  of  ARC  and  NODE,  respectively,  to 
test  the  hypotheses  in  Table  1. 


Computation  of  Covariance  Matrix 

The  estimation  of  the  covariance  matrix  under  HQ,  H-j,  H,,  and  U3  uses 

(12)  e 9 L 8 In  L _ n 3f(x,a,b,c)  , 8 f(x,a,b,c)  , 1_ 

8 8 8 <f>  j '0  30  8$  f(x,a,b,c) 


together  with  the  expressions  in  Table  1 In  ESTIMA  and  HYP123.  These  sub- 
routines employ  numerical  integration,  as  described  in  [12,  p.923]  to 
evaluate  l , using  a,  b,  c in  place  of  a,  b and  c respectively.  Although 
the  weights  in  the  W and  Y arrays  apply  for  double  precision  computation, 
experience  has  shown  little  loss  of  accuracy  by  using  single  precision. 

Figure  2 offers  a>.  example  of  the  output  for  100  observations  drawn  from 
f (x,l ,5,12) . 


. ■ * < 
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Likelihood  Ratio  Test 

Since  parsimony  clearly  has  advantages  in  modeling,  one  wants  to  test 
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the hypotheses  in  Table  1 to  see  if  one  or  two  parameters  can  be  elimi- 

A A A 

nated  from  the  representation  (1).  Let  LQ(,  a^»  b..,  c^)  denote  the 
maximum  of  the  likelihood  function  under  K..  where  i = 0 corresponds 
to  (1)  and  X = (Xj,...,Xn).  For  example,  LU.a-j »b-j  ,c-| ) = L(X,0,b-j  ,c-j) 
and  L{X.,a2,b2,c2)  = L(X,a2>a2,c2).  Then  the  likelihood  ratio 


(13) 


N A A (AAA 

Rj  = L(X^  ,b.j ,c^ )/L(X.,<JQ»b^,CQ)  i - !,*••> 6 


lies  in  (0,1).  The  closer  R-  is  to  rnity  the  more  credible  is  the 
hypothesis.  Although  the  distribution  of  under  is  beyond  our 
reach  it  is  know  the-  as  n increases  the  distribution  of  -2  £n  R^ 
converges  to  the  chi-square  distribution  with  degrees  of  freedom  equal  to 
the  number  of  constraints  imposed  by  the  hypothesis  [10],  For  H^,  Hg  and 
there  is  1 degree  of  freedom;  for  and  there  are  2 degrees 

of  freedom.  Therefore 

-XfO-«)/2 

(14)  pr(R,j  * e ) = 1 - a 


where  XfO-a)  denotes  the  1 - a critical  value  for  f degrees  of 
freedom  and  f a 1 for  i 3 1... .,3  and  f = 2 for  1 4,..., 6. 

Table  2 shows  critical  values  of  R^  corresponding  to  tests  of  selected 
siies. 
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MAXI  MUM  LIKELIHOOD  PSTIH ATION 

P(X)  =G(X,  A,B,C)  fG  (X,  D,A»C)+G  (X#C,A#B) 

G<X,T,P,R)  =T*EXP  (-X/T)  /(  (T-P)  *{T-R) ) 

SAMPLE  MEAN*  0.209 147E  02  SAHPL?  ¥ABIANCF=  0.250988F  03 
DELTA*  0.697155E-02 

A-  0.474065F  00  B=  0.S5S364E  01  C*  0.148639S  02 


-r&H 


0.  137369B  01 


0.  100000*  01 


COVAPIA  PCF  MATRIX 
-0.359012E  01 
0.  195890S  02 

COPRKtATION  MATRIX 
- 0. 692083*  00 
0.  100000B  01 


0.30U326E  01 
•0.  22248 IE  02 
0.  316986E  02 

0.461184P  00 
■0. 8Q2826P  00 
0.  100000E  01 


HYPOTHESIS  1:  A«0,  B<=C 


I*;.1 


1 

'.Vi 


8=  0.646436E  01  C-  0.  144503K  02 

VA8<a)  o 0.17Q784B  02  VAR  <C>  = 0.  36  14  12P  02  COV  (B.C)  »*0.231085E  02 

COBB  <8#C)»-0. 9 30136P  00 
LIKELIHOOD  PATIO*  Q.781658P  00 


HYPOTHESIS  2:  A»8<»C 


8-  Q.256668F  01  C*  0.1S7813P  02 

V AS  < PI  3 0.  $75024*  00  VAP(C)«  C.  1029190  02  COV (B,C| »-0. 171543B  01 

COP9(R,C)*-0.705151P  00 
LIKELIHOOD  RATIO-  0.72256 3P  00 
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HYPOTHPSIS  3: 


A<*  S-C 


A*  0.27 2 369K-0 1 C-  0.104437F  02 

V?  R{S)  = 0.I55362F  00  VAR (C) = 0.9508158  00  CO?  (A,C)  = -0.  172593F.  00 

COPP(A,C)=-0.349580F  00 
LIKELIHOOD  RATIO*  0.575515F  00 


HYPOTHESIS  4:  A=B=0 

C=  0.209147P  02 

.95  LONEP  POINT”  0.172517.  02  .95  UPPER  POINT*  0.257062E  02 

LIKELIHOOD  RATIO*  0.198Q47F-04 


HYPOTHESIS  5:  A=0,  B=C 

C=  0.1Q4573R  C2 

.95  LOWEP  POINT-  0.91466SF  01  .95  UPPP~  POIKT*  0.120730F  02 

LIKELIHOOD  RATIO*  0.5657638  Oo 


HYPOTHESIS  6:  A = B=C 

C=  0.697155?  Cl 

.95  LOWER  POINT*  0. 524519?  01  .95  UPPEP  POUT*  0.783313P  01 

LIKELIHOOD  RATIO*  0.67690 1F-0 3 


AVI 


Figure  2 (continued) 
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Table  2 

Critical  Values  of  R^  for  Tests  of  Selected  Sizes 


In  the  case  of  2 d.f.  the  chi  square  distribution  is  E{1).  There- 
2 

fore,  R.j  is  the  probability  that  under  H,j  (ia4,5,6)  one  would  observe 
a likelihood  ratio  less  than+  R^.  For  example,  under  Hg  in  Figure  2 
R5  = 0.5658  and  r|  = 0.3201 . 

Confidence  Intervals 

Let  us  first  concentrate  on  Hg  and  Hg.  Under  c^  = X/(i-3) 

A 

and  n c^/c  has  the  chi-square  distribution  with  (i-3)n/2  degrees  of 
freedom.  The  ERLANG  subroutine  uses  this  fact  to  compute  a confidence 
interval  for  c and  relies  on  the  CHISQ  subroutine  to  provide  critical 
values  of  chi-square. 

For  Hq,  Hj,  H2  and  H3  no  similar  theory  is  available.  However,  if 
n is  sufficiently  large,  one  can  compute  approximate  individual  confidence 
intervals  for  a,  b and  c,  using  the  estimated  variances  in  the  corresponding 
covariance  matrix.  Experience  with  the  set  of  subprograms  in  the  appendix 
has  revealed  that  even  for  n - 100  the  sample  var(a),  var(b),  var(c)  are 

A A A 

large  relative  to  a,  b and  c respectively. 


+Here  R^  is  called  the  P-value.  See  [8]. 
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Bias  Considerations 

A A A 

In  small  and  moderate  size  samples  a,  b and  c are  biased.  In 

A A 

particular,  experience  has  shown  that  a overestimates  a and  c 

A 

underestimates  c.  Since  c does  most  to  affect  the  shape  of  the  tail 
of  the  distribution  we  especially  want  to  consider  ways  of  reducing  bias 
for  this  quantity.  One  approach  to  bias  reduction  uses  the  jaokknife 
method  [11]. 

The  elementary  form  of  the  jackknife  method  removes  bias  to  order 
1/n.  Suppose  c is  computed  using  n observations  and  c^  and  c^ 
computed  using  the  first  m = n/2  observations  and  the  last  in  = n/2 
observations  respectively.  Then  one  can  easily  show  that 


(15) 


c»  2c  - (cO)  + c^)/2 


is  free  from  bias  to  order  1/n.  Notice  that  the  computation  of  c 
requires  3 passes  through  the  estimation  procedure. 

More  powerful  jackknife  methods  of  bias  reduction  are  available  [11] 
Our  reluctance  to  incorporate  any  one  of  them  into  the  estimation  pro- 
cedure is  a consequence  of  the  additional  cost  they  imply.  However,  a 
user  of  the  estimation  procedure  in  the  Appendix  can  easily  write  a bias 
reduction  program  to  use  in  conjunction  with  ESTINA. 


* 
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4.  Appendix 


SUBROUTINE  ESTIMA(X,N,NUH) 

C 

CO  HUE NT  THIS  SUBROUTINE  CONDUCTS  A GRID  SEARCH  ON  A IN  INCREMENTS  OF 

C DBLTA 


INTBGEP  I,J,K,N,N,NUM 

REAL  A.B,C,AS,BS,CS,AA(6),BB(6)  ,CC(6)  ,CORR(3,3)  ,COV  (3,3)  ,D(3,3|  , 

2 DELTA,  DEN,  F,H  (3)  , LC,LIKE  (6)  ,LLP,LOGX,HAXLLP,UC,tl  (15)  ,X(N)  , 

3 XBAR,XSUN,Y(15) 

DATA  »/. 2345781703,, 560 1008 428,. 9870082629, 1. 22 3664 40 215, 

2 1.57444872163,1.94475197653,2.34150205664,2.77404192683, 

3 3.25564334640,3.80631171423,4.45847775384,5.27001778443, 

4 6. 35956346973,8.03178763212, 11.5277721009/ 

DATA  X/. 093 3078 120,. 49 26 91 7403, 1.2155954121,2.2699495262, 

2 3.6676227218,5.4  253366  27  4,7.5659162266,10.1202285680, 

3 13.1302824822,16.6544077083,20.7764788994,25.6238942268, 

4 31.4075191698,38.5306833065,48.0260855727/ 

1 PORHAT(1H1,25X»MAXIHUM  LIKELIHOOD  ESTIMATION* 

2 *//22X,  • P (X)  -G  (X , A, B, C)  fG(X,  B,A,C)fG  (X,C,A,B)  »//22X 

3,*G(X,T,P,P)=T*EXP(-X/T)/((T*P)*(T-P))»//4X,*N=«,I5,«  SAMPLE  M 
4BAN=*,E13.6»*  SAMPLE  VARIANCE-* , E13.6//30X, *DELTA=* ,E13.6//12X 
5,  *A=*,B13.6,*  B»*  ,P1 3*  6 , ' C = * ,E13.6///) 

2 FORMAT  (•  SPE  HYPOTHESIS* , 12/////) 

3 FORMAT  (3 IX , *COVAPIANCP  MATRIX* //15X,3 (El 3. 6,5X) //3 3X, 2 (E13.6,5X) / 
2/51 X,E 13. 6//31X, ‘CORRELATION  MATBIX*//1SX,3  (B13.6, 5X) //33X, 2 (E13.6 
3 , 5X ) //5 1 X , E 1 3.  6/////) 

XSUM*0 

LOGX=0 

DO  100  I*  1, N 
LOGI=LOGX+ ALOG (X (I) ) 

100  XSUM**XSUN’fX  (I) 

XBAP*XSUM/N 

A=0. 

DBLTA^XBAP/  (3.  *NUH) 

M-NUH-1 

LLP*0 

DO  150  T*  1 , M 

HAXLLP»LLP 

AS=*A 

BS=B 

cs=»c 


t 

This  set  of  FORTRAN  subroutines  computes  the  maximum  likelihood  estimates  of 
a,  b and  c in  f(x,a,b,c)  for  Hq  through  Hg  in  Table  1.  X denotes  the 
floating  point  data  array,  N denotes  the  sample  size  and  NUN  denotes  the 
resolution  DELTA  * X/(3*NUN)  for  conducting  the  grid  search. 


»r«rsy.- 


-18- 


APC  AND  NODE  SEARCH  ON  THE  BOUNDARIES  OF  THE  FEASIBLE  FFGION 


LC=  (XBAR-A)  /2. 

UC=XBAR-2.* A 

CALL  UPDATF  (X, N , XBAF ,LC ,UC,  A , B,  C,DSLT A,  1,LLF) 

IF  (SAXLLF. FQ. 0.)  M AXLL  F=LLF 
IF  (BAXLLF. GT.LLF)  GO  TO  160 
150  A=AfDELT A 
160  A=AS 
B=BS 
C=CS 
C 

CONSENT 
C 

DO  170  1=1, 3 

CALL  AFC  (X,N,XBAP,A,B,C,DBLTA,MAXLLF,I,AA(I),BB(I>,CC(I)  .LIKE  (I) ) 
170  CALL  NODE(N,XBAB,LOGX,A,B,C,HAXLLF,I, 

2AA  (I  + J)  ,BB(I+3)  ,CC  <1-4-3 ) ,LIKE(U3)) 

C 

CONSENT  OUTPUT  COMPUTATIONS  FOLLOW 
C 

SSQ=0 

DO  180  1=1,3 
DO  180  J = I,  3 
180  D(I, J) =0 

DO  190  T=1,N 

190  3SU=SSQ4-  (X  (I)  -X  BAT  ) **2 
SSQ=SSQ/(N~1) 

WRITS  (3,1)  N ,X  BAF , SSQ ,D2LTA,A,B,C 

1=0 

IP  (A. FQ.  0. AND. B. LT.C)  1=1 
IF  (A. BO. B. AND. B. LT.C)  1=2 
IF  (A. LT.  B.  AND.  B.EO.C)  1 = 3 
IF  (A.EQ.O. AND.B.FQ.O)  1=4 
IF  (A. GO.O, ANT. B.EO.C)  1=5 
IF  (A.FQ.B. AND. B.EO.C)  1=6 
IP  (I.FO.O)  GO  TO  200 
WPITE  (3,2)  I 
GO  TO  450 

200  DO  300  1=1,15 

CALL  COHPUT (X (I) , 1 , B,C , A , 1 , H (1)  ,F) 

CALL  CCNPUT  (X  (I ) , 1 , A,C, B , 1,  H { 2)  , F) 

CALL  COM  PUT  (X  (I)  ,1,A,B,C,1,H(3)  ,P) 

F=EXP  (F) 

DO  300  J=1, 3 
DO  300  K= J , 3 

300  D(J,  K)  =D  («],  K)  4-H  (J)  *H  (K)  *F*W  (I) 

DEN=  (D  (1 , 1)  *D  (2  ,2)  *D  (3, 3)f2 . *D  ( 1, 2)  *0(2, 3)  * D ( 1,  3) 

2-0(2, 2)  *D(1,3)**2-n(3,))*D<1,2)  **2-P<1,1)*D  (2,3)**2)*N 


«^K5£2a?SSwe^^  afe^M. 
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C0V(1,  1)  = (D  (2,2)  *D  (3,3)  -D  (2,3)  **2)  /DEN 
COV(1,2)  =(-D(1,2)*D(3,3)+D(1,  3)  *D(2,3))  /DEN 
COV(1,3)  = (D  (1,2)  *D  (2,3)  -D  (1 ,3)  *D  (2,2) ) /DEN 
COV  (2, 2)  3 (D  (1,  1)  *P(3,3)  -D(1,3)**2)/DEN 
COV  (2,3)  = (-D(1,1)  *D  (2,3)  4-D  { 1, 2)  *D(1,3))  /DEN 
COV  (3,  3)  = (D  (1,  1)  *D  (2,2)  -D  (1 , 2)  **2)  /DEN 
DO  400  1*1,3 
DO  400  J*I,3 

400  CORR  (I,  J)  =COV  (I,  J)  /SQRT  (COV  (1,1)  *C0V  (J,  J) ) 

WRITE  (3,3)  COV(1,1)  ,COV  (1 ,2)  ,COV  (1 , 3)  ,COV  (2 , 2)  , COV  (2,  3)  , 

2COV  (3,3)  , CORR  (1,1)  ,CORR  (1,2)  ,CORR(1, 3)  ,COPR  (2,2)  ,COFR  (2,3)  , 
3COFR  (3,3) 

C 

COMMENT  CHECK  HYPOTHESES  1,2  AND  3 
C 

CALL  HTP123 (X, N,XBAR, LOGX, DELTA ,H  AXLLP,  AA,B B,CC , LIKE) 

C 

COHMENT  CHECK  HYPOTHESES  4,5  AND  6 
C 

DO  500  1=1,3 

500  CALL  ERLANG  (N,XBAR, LOGX, HAXLLF,I) 

END 


DOUBLE  PRECISION  FUNCTION  G (Y, THETA, PHI, RHO) 

C 

COHMENT  COMPUTES  THETA*BXP  (-Y/THBTA  )/(  (THETA-PHI)  * (THET  A-RHO) ) 
C 

PEAL  PHI, RHO, S, THETA, Y 
REAL*8  ARG, CHECK, Z,ZZ,ZZZ 
G=0 

IF  (THETA. EQ.O.)  RETURN 

ARG=Y/THETA 

S=1. 

Z=THBTA/$  (THETA-PHI)*  (THETA-RHO) ) 

IF  (Z.LT.O.  ) S3-S 
IP  (ARG. LE.  174. 673)  GO  TO  25 
10  ARG=-ARG+DLOG (DABS  (Z) ) 

IF  (A8G.LT.-1D0.218)  PETUBN 
G=S*D8XP  (ARG) 

RETURN 

25  ZZZ*DEXP  (ARG) 

CHECK3 (10D-78) *ZZZ 
ZZ-DABS  (Z) 

IP  (ZZ. IT. CHECK)  RETURN 

G=Z/ZZZ 

RETURN 

END 


t.&«i  • s^- 
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S0BROUTINE  UPDATE (X,N, XBAR,LC,UC, A, B,C„ BELT  A, J, LLP) 


COHHENT 

C 


PERFORMS  BINARY  SEARCH  ON  C FOR  GIVEN  A 


INTEGER  J , N 

REAL  A,B,C,CC, DEL, DELTA , DERI VC, LC,LLF,HC,H(12),X(N)  ,XBAR,V(16| 
DATA  V/0.,0.,0.,1.,1.,0. ,0. ,0.,0.,0.,1. ,0., 0.,0.,0. ,-2./ 

DATA  H/1. ,1 .,.5,0.,-1.,0.,0.,0. ,-1. 

CC=  (LC-fOC)  /2» 

100  c=cc 

B=XBAR*W  (J)  +A*H  (J+4)  +C*W  (04-8) 

A=XBAR*V(J)  +A*VCJ+4HB*V  (j+B)  .|.C*V{J+12) 

CALL  CONPOT  (X,N, A, B ,C, J, DBRIYC, LLP) 

IP  (DBRIVC.GB.O)  LC=C 
IP  (DBRIVC j LE.  0)  OC=C 
CC=  (LC-fUC)  /2« 

DEL  = ABS  (C-CC) 

IP  ( DEL. GT- DELTA)  GO  TO  100 

RETURN 

END 


SUBROUTINE  COHPUT(Y,M,T,P,R,J,DEPIVC,LLr) 

C 

COMMENT  COMPUTES  LOGLIKBLI HOOD  DERIVATIVE  WITH  RESPECT  TO  C 
C 

INTEGER 

REAL  DBPTVC,LLP,P,R,?,Y {B) 

RBAL*8  P,G,GP,GP,GT 

LLP-0 

DEBIVC=0 

GO  TO  (100,200,300,400),  J 
100  DO  150  K=1 , N 

GTaG  (X  (K) ,T,P,F) 

GP=G  (I  (K)  ,P ,T,R) 

GR=G  <Y  (K) , P,T,P) 

P=GTK»P4-GR 

LLF*LLP+DLOG(P) 

GT*GT/P 

GP=GP/P 

GR=t5R/P 

150  DBPXVC=DEPIVCK»P*(1./PPV  (K)  /R**  2+1. / (P-NH-1  ./<"*-P)) 

2+GP/  (P-F)-HST/(T-R) 

RETURN 

200  DO  250  K*  1 , H 

GP*G  |T  (K)  ,P,R,0) 


swag- — ^ . 
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gr=g  a (ki  ,p,p,oj 

F=GP+GP 

LLF=LLFfDLOG  (P) 

GP~G  P/F 

gp=c;r/f 

2 50  DBF  I VC=PEPI  VC+GP/ (P-F)  *GP*  { Y (K)  /R**  2- 1 . / <R- P) ) 

RETURN 

300  DO  350  K-  1 , M 

GP=G  (Y  (K > ,P,P.,F) 

GB»G  «Y  f K)  ,R,P,P) 

F=GFPGP*  ( (P-F)  $ Y (K) -P*»)  /P**2 
LLF  = LLP+BLOG(F) 

GP=GP/F 

GF=GR/P 

350  DBF IVC= D?PT VC+GF * (1  ,/F  + Y (KJ  /?  *♦2-2.  / (B-P) ) 

2fGP*  ( (P~P)  *Y(K)+P*  (P4-P)  > /(  (R-P)  *P**2) 

, RETURN 
400  DO  450  K=1#N 

G?=G  (Y  (K) ,T.F,F) 

GE=:;  JY(K>  ,R,T,T) 

F=6TH»P*  ( (f-T)  *Y  (K| -T*F  ) /p**2 
LIF=  LLFf DLOG (?) 

GT  aGT/F 
GRaGR/^ 

4 50  OEF I VC  a DEF I VC  f 2 . *G?/ (T- P ) 

2 4-G P « ( ( (P*T)  *Y(K)  *T*P)  * ( Y (K)  /E*«2-1.  /R-2,/'F~T) ) +Y(K)-T)  /S**2 
RETURN 
END 


C 


SUBROUTINE 


ARC(X,N,XBAP,A,H,C,0<5LTA,KAXLLF,  I,  AA,  BB,  CC.LIKF) 


CQDNENT  COMPOTES  SOLUTIONS  POP  APCS  AND  APPLIES  TO  HYPOTHESES  1-2  AND  3 
C 

INTEGER  I » N 

BEAL  A,  AA,P,P8»CfCC,R^LTA,I  INF, H AXI  t?,U (6) , I (N) #X8*R 
DATA  U/.  5» . 

CALI  UPDATE  <X,N#XPA*.XBAP*tt(X»  0 XBAR*U  (If  3|  # AA» BB»CC, DELTA  , If  1 , 
2LIKE) 

H (NAXLIP.G^.LIKE)  »FT'IPN 

A?  A A 

U-OU 


C^CC 

KAXLLF^LIKP 

RPTURS 

END 


u u u 
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SUBROUTlNF  NODF  (N,XBAF  f LOGX ,A ,B ,C ,M AXLLF, I, AA , B B,CC ,LIKF) 

0 MMENT  COMPUTES  SOLUTIONS  FOP  NODES  AND  APPLIES  TO  HYPOTHESES  4,5,6 
INTEGER  I ,N 

PEAL  A, AA,B,BB,C,CC,LTKE,lOG2 ,LOGX, NAXLLF,XBAR,  V (9) 

DATA  W/0. ,0.  ,.  3333  33,0.  ,.5,  .3  33  333,  1. , . 5, . 3 3333 3/, LOG2/. 693 147/ 
LIKE=-N*  (I*  (1.  + ALOG  (XBAF/I)  )+3.  *H  (I)  *LQG2)  + (1-1)  *LOGX 
AA=XBAR*W(I) 

BB=XBAF*W  (1+3) 

CC  = XBAF*  W (I  +6) 

IF  (LLF.GT.LIKE)  FFTUFN 

A = AA 

B=BB 

C=CC 

MAXLLF=LIKF 

FKTUPN 

END 


SUBROUTINE  H YP 1 2 3 ( X, N, XB AP , LOGX , DELTA.MAXLL P, AA, BB, CC,LIKF) 


COMMENT 

C 


PERFORMS  OUTPUT  ANALYSES  POP  HYPOTHESES  1,2  AND  3 


1 

2 

3 

4 


INTEGRA  I , J ,K  (°) ,KA,KB,KC,L,N 

REAL  A,  A A (6),8.BB(6)  ,C,  CC  (6)  , CP  B,CCC,CBC,  CC  BB,CCCC,  CCBC  , DFLTA  , 
0(2,2)  , DPN , F, H B, HC, LIKE (6)  ,LOGX,LPATIO ,MAX LLP, U (6) , W ( 18) , 

X (N)  , X B A P , Y (15) 

DATA  K/1 * 1,2, 2, 3, 3 ,1,2,1/ 

DATA  U/.  5, .3  33  3 33, .3313  13,1.,  1.  , . 5/ 

DATA  «/.  239  5?81703,.5601C08428,  .8870082629,  1.22  3664  40  215, 

I. 57  4448  7216  3,1.94475197  653,2.  JU  1502 0566  4, 2.  77404  19268 3, 

3.  25 664  1)4640,  1. 8C631171 42 3, 4. 45847775384,6.  2700 177846  1, 

6.  3c95634697.),e.03lf878J212,  1 1. S27-*?21009/ 

DATA  Y/.0931078120,. 4926917403,  1. 21 5598412 1 ,2.2  619495  26  2, 

).f  6 76227218,5.4  25336627  4,7.  5659  162266,  1 0. 1202285690, 

II. 1  1029  2 !» 82 2, 16. 65440 77 00 3, 20. 7 7647 8099 4, 25. 6239 9 42 26 8, 
n.4  07M<M6«fl,  18, 5 10  693  306  5, 46. 0 2f  08  55  727/ 

FORMAT  {•  HYPOTHESIS  1:  A*0,  B<=C*//1 

FORMAT  (•  HYPOTHESIS  2:  A = F<=C'//) 

FORMAT  (•  DIYPOTHESIS  3:  A<  = 3=C  •//) 

PORH\T(15X,  »A=».E1  J.  6,8  X,  * B = * , E 13. 6 , 5X,  *C=*  ,E13.6//*  SFP  HYP 

20THKS IS  ' , I 2/////) 


2 

3 


5 


...  », 
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FOB  NAT(22X,  • B=*  ,E13.6,5X,*C=*  ,E  1 3.6//9X, « VP.  B (B)  = * ,F  1 3.  6, 5X,  • VAP  (C) 
2=',E13.6,5X,*COV(B,C)  =• , El  3. 6//28X, *COPR <B, C) -• , B 13. 6//25X, 'LTKELI 
3H00D  PATIO= ' , El  3. 6/////) 

6 PORWAT  (22X, «A=* ,S1 3.6,5 X, *C=* ,F 13. 6//9X, 'VAR (X) = • , E 1 3.6, 5X. » VAP (C) 
2=«,E13.  6,5X,«COV (A,C) =• ,E13.6//28X, *COPR<A,C)  = » ,B13.6//25X, *LIKELI 
3HOOD  RATIO= • ,B1 3.6/////) 

DO  500  1=1,3 

IF  (I. BQk 1)  WRITS  (3,1) 

xp  (I.  BO.  2)  WRITS  (3,2) 

IP  (I.BQ.3)  WRITE  (3,3) 

L=0 

KA=K  (I) 

KB  = K (14-3) 

KC=K  (14-6) 

L=I 

DO  100  J=KA,KB,KC 
JJ=Jf3 

IP  (LIKE  (I)  .IT. LIKE (JJ)  ) L=JJ 
A=AA  (L) 

B = BB  (L) 

C=CC  (L) 

100  LIKB(I)  = LIKB(L) 

L=0 

IF  (A.BQ.O.  AND.B.EQ.O)  L=4 
IP  (A.BQ.O.  AND.  B.BQ.C)  L=5 
IF  (A.  BQ.  B.  AND.8.EQ.C)  1=6 
IF  (L.LT.4)  GO  TO  150 
WRITS  (3,4)  A,B,C,L 
GO  TO  500 

150  LFATIO=iiXP  (l IBB  (I)  -KAXILP) 

DO  175  J= 1, 2 
DO  175  L=J, 2 
175  D (J,L)  =0 

DO  475  J=1,  15 
GO  TO  (200,300,400),  I 
200  CALL  CONPUf(Y(J|  ,1,  A,C,  B,  2,  HB,  F) 

CALL  C0HPUT(Y(J),1,A,B,C,2,HC,F> 

GO  TO  450 

300  CALL  CQHPttT  (Y  (J)  ,1,C,B,A,4,HB,F) 

CALL  COHPUT  (X  (J)  ,1 , l,B,C,  3,  KC,P) 

GO  TO  450 

400  CALL  COHPOT  (Y  (J) ,1,C,B,A*3,H9,F) 

CALL  COHPOT  (Y(3),1,A,B,C,4,KC,P) 

450  P=EXP (F) 

D(1,  1)  =D(f,  1)+HB»*2*r*W  (J) 

D (2,  2)  =D  (2, 2)  +HC**2*F*W  (J) 

475  D(1, 2)  =D  (1, 2)  4H8*HC*P*W  (J) 

DEN  = (D  (1 , 1)  *D  (2 , 2)  -D  ( 1 , 2)  **2)  •« 


•C-  W t'J 


-24" 


CBB=0(2,2)/DSN 
CCC-D  (1,1)  /DEN 
CBC=-D  (1,2)  /DEN 

CCBC=-D(1,2)/SQRT(D(1,1)  *D(2,2)  ) 

IF  (t.EQ.  1)  WRITE  {3,5)  B,C,CBB  ,CCC  ,CBC,CCBC,LRATIO 
IF  (X.EQ.2)  WRITE  (3,5)  B,C,CBB,CCC,CBC,CCBC,LRATIO 
IF  (X.EQ.3)  WRITE  (3,6)  A,C  ,CBB  ,CCC  ,CBC  ,CCBC,LP  ATIO 
5C0  CONTINUE 
RETURN 
END 


SUBROUTINE  ERLANG  (N  ,X BA  » ,LOGX  , I*  AXLI.. F, I) 

C 

COMMENT  PERFORMS  OUTPUT  ANALYSES  FOR  HYPOTHESES  4,5  AND  6 
C 

INTEGER  I ,N 

SEAL  C, CHIS Q, DF, LC, LLF,LOG2 ,LOGX,LR ATIO, MAX LLF, UC,P  (3),XBAR 
DATA  W/0.,0.,1./,  LOG2/. 693147/ 

1 FORMAT  (°  HYPOTHESIS  4:  A=B=0'//) 

FORMAT  (•  HYPOTHESIS  5:  A=0,  B=C»//) 

FORMAT  (•  HYPOTHESIS  6:  A=B=C‘//) 

FORMAT (3?X, *C* * ,E13. 6//8X, * ,9 5 LOWER  POIKT= *,E1 3. 6, 5X, • , 95  UPPER  P 
20INT^=  * , & 1 3.  6//25X,  • LIKELIHOOD  R ATIO=*  , E 1 3. 6 /////) 

C=XBAR/I 
DF~2. *I*N 

LC=DF*C/CHISQ  (DF,.975) 

UC=DP*C/CHISQ  (DP, .025) 

LLF=-N*  (I*  ( 1.  mOG  (C) ) 4-W  (I)  *LOG2)4-  (I- 1)  *LOGX 
LRATIO=EXP (LLP-HAXLLF) 

IF  (I.  EQ*  1)  WRITE  (3,1) 

IF  (I.EQ.2)  WRITE  (3,2) 

IF  (I.EQ.3)  WRITE  (3,3) 

WRITE  (3,4)  C,LC,UC,LRATI0 
RETURN 
END 


FUNCTION  CHISQ(DP,P) 


CONSENT  COMPUTES  CRITICAL  TALUS  OF  CHI-SQUAF2  FOP  PPOBA BILITY  P 
C AMD  DF  D2GBEBS  OF  PHEEDOH 

INTEGBH  I 
REAL  C (3)  PD  «3| 

RB1L  DF,P,Q,T, XP,l?0?!,DBN,f  ,SQDF  ,SQHALF,  Y2,Y3,Y4  PY5,  Y6,Y7.  H(7) 
DATA  C/2. 515517,. 802853,. Q10328/,D/1. 432788, .189269. .001308/ 
HUH*0 
DEH=1. 

Q*P 

IF  (P.LE..5)  GO  TO  5 
Q=1.-P 

5 T=SQBT  (ALOG  (1./Q**2) } 

DO  10  1=1,3 
HDH  = IIUM+C  (I) 

10  DBN=DBN+D(I)*T**I 
XP-T-MUB/DBN 
IF  (P.GB..S)  GO  TO  15 
XP=~XP 
15  Y=XP 

SQDF=SQRT  (DF) 

SQ8ALF»SQPT(.5) 

Y2»Y*Y 

Y3«Y*Y2 

Y4=Y3*Y 

Y5*Y4*Y 

Y6=Y5*Y 

Y7=Y6*Y 

H f 1|  3Y/$QHALP 

H (2)  = 2.*<Y2-1.)  /3, 

tt(3J  = (Y3-7.  *Y)*SQHALF/9» 

H (I)  *-  f6*Y4+1A.  *Y2-32.)  /405. 

H(5)  =(9.*Y54-256.*Y3-t»33.  *Y)  *SQH ALP/48S0. 

H (6)  = ( 12. *Y6-24  3„  *vu-923.*Y2f  1472.)  /25515. 

H<7)  =-  (3753. *Y74-43M.*Y5- 289517. •yi-289717.  *Y)  •SQHALF/9105400. 
CHISQ=1. 

CO  20  1=1,7 

20  CHISQ=CHIS04-H(I)/S0DP**I 
CHISO*CRISO*DF 
RETURN 
BHD 
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